Numerical study of drying process and columnar fracture process in granules-water 

mixtures 

Akihiro NishimotcQ 

Department of Physics, Kyoto University, Kyoto 606-8502, Japan 

Tsuyoshi Mizuguchi 

Department of Mathematical Sciences, Osaka Prefecture University, Sakai 599-8531, Japan 

So Kitsunezaki 

Graduate School of Human Culture, Nara Women's University, Nara 630-8506, Japan 

(Dated: February 6, 2008) 

The formation of three-dimensional prismatic cracks in the drying process of starch- water mixtures 
is investigated numerically. We assume that the mixture is an elastic porous medium which possesses 
a stress field and a water content field. The evolution of both fields are represented by a spring 
network and a phenomenological model with the water potential, respectively. We find that the 
water content distribution has a propagating front which is not explained by a simple diffusion 
process. The prismatic structure of cracks driven by the water content field is observed. The 
depth dependence and the coarsening process of the columnar structure are also studied. The 
particle diameter dependence of the scale of the columns and the effect of the crack networks on the 
dynamics of the water content field are also discussed. 

PACS numbers: 62.20.Mk, 46.50.+a, 45.70.Qj, 47.56. +r 



I. INTRODUCTION 

Crack patterns are observed in everyday life 
Columnar joint in cooled lava, especially, is an intriguing 
phenomenon and has fascinated many people for cen- 
turies and has been studied by field works 0, d, 0, HI- 
Theoretical studies have been carried out from the view- 
point of the analysis of the thermal field @, 0, @|, the 
ordering of patterns 0, EH , and size selection of columns 
by the finite element method (Tlj . However, many ques- 
tions remain open, such as, 'why are the prismatic struc- 
tures with the polygonal cross sections formed?', 'how 
are the size of columns and the distribution of polygon 
types determined?'. The scale of columns is 0.1 ~ 1 m 
in diameter, hence it is too hard to study the fracture 
process of cooled lava experimentally. 

Recently, the patterns of cracks similar to columnar 
joint are investi gate d in the dryi ng p rocess of starch- 
water mixtures [H, El El El M, El E3- Well- 
controlled experiments are possible due to the small scale 
of columns (1 ~ 10 mm in diameter). In typical experi- 
ments of these studies, after the mixtures of starch and 
water are poured into a container, water evaporates from 
the sample surface. During the drying process there are 
three stages [H EE El EH- Stage 1: The water con- 
tent decreases uniformly. Stage 2: Cracks are formed 
to extend to the bottom of the sample and a sudden 
change of water content occurs. These cracks have a 
uniform structure along the vertical direction and are re- 
ferred as the primary (type I) cracks [jj, El El, HI H3 • 
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In the case that the thickness of a sample is sufficiently 
large, the primary cracks appear only between the mix- 
ture and the side wall of the container. Stage 3: The 
water content decreases slowly and the water content 
distribution becomes non- uniform, i.e., it has a "front" 
[T7L E|[ . The front propagates inward and the mixture 
shrinks non-uniformly and another type of cracks, i.e., 
the secondary (type II) cracks are formed. The secondary 
cracks show three-dimensional prismatic structure. The 
size of columns depends on the evaporation rate [H EH 
and it coarsens with the distance from the sample surface 
[HE3- I n this paper, we focus on the secondary cracks. 

By a numerical model, Hayakawa [13] exhibited that a 
regular prismatic crack structure is formed by the steady 
sweep of external field (temperature) with fixed shape. 
Experimental results of the drying starch, however, sug- 
gest that the shape of the external field (water content) 
changed temporally. In order to study the crack pat- 
tern of starch-water mixtures, the drying process and 
the fracture process should be taken into account three- 
dimensionally. Combined numerical analysis of both pro- 
cesses, however, has not been performed yet. The drying 
process involving all stages, especially, has not been an- 
alyzed, that is mainly because the water transportation 
in the porous media are complicated process. 

The non-uniform shrinkage due to the change of the 
external field causes directional propagation of cracks, 
and the cracks become the boundary condition of the ex- 
ternal field, simultaneously. The changes of the external 
field can be induced by cracks as the boundary condi- 
tion, as suggested by some studies [H [H |25I |. This 
type of crack is referred as self-driven crack. It has nei- 
ther been discussed in detail nor demonstrated numeri- 
cally whether this effect is important or not in the case 
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FIG. 1: Schematic diagram of the model. Three elements 
and interactions between them are drawn. The interactions 
indicated by solid arrows are taken into account in our model 
while the dashed ones are ignored. 



of starch-water mixtures. In some studies on columnar 
joints 0,0, SHE S3 and desiccation cracks [13, SI, the 
effect of cracks as the boundary condition is considered 
to play an important role. On the other hand, in studies 
such as two-dimensional directional fracture of glass strip 
29|. l30l . l3llh c olumnar joints [1] and desiccation cracks 
12L fla. [l8ll32l. [33l [34l|. this effect is hypothesized not to 
be important. 

In this paper, our goal is a comprehensive descrip- 
tion of the water content distribution and the three- 
dimensional realization of the crack formation process 
driven by the non-uniform shrinkage during the drying 
process of starch-water mixtures. 

This paper is organized as follows: in Sec. II, details 
of the model are described. We assume that the mix- 
ture is an elastic porous medium which has a stress field 
and a water content field. The evolution of the former is 
represented by a spring network model and the latter is 
represented by a phenomenological model with the water 
potential. Here we suppose that cracks do not affect the 
evolution of the water content distribution. In Sec. Ill, 
numerical results are presented. A propagating front of 
water content distribution and the prismatic pattern of 
cracks driven by the water content field are obtained, 
which well capture those observed in the experiments. 
The particle diameter dependence of the scale of columns 
is suggested. In Sec. IV, the properties of the drying front 
and the validity of some assumptions used in our model 
are discussed. The effect of the crack networks on the 
dynamics of the water content field and the interaction 
between the stress field and the water content field are 
considered. The geological columnar joints and the char- 
acterization of the crack patterns are discussed. Finally, 
we conclude this paper with a summary in Sec. V. 



II. MODEL 

The drying fracture process consists of various factors. 
Here, as key factors, we focus on the following quantities 
and processes (Fig. [1} : the change of water content by 
the evaporation and transportation of liquid water and 
vapor in porous media, the non-uniform volumetric con- 
traction and the stress intensification, the formation of 
cracks, and the change of the boundary condition of the 
stress field. 

Starch-water mixtures are elastic porous media com- 
posed of solid phase (granules), liquid phase (water) and 
gas phase (air and vapor). As the water content decreases 
by drying, the state of the mixture changes from "capil- 
lary fringe" (liquid is saturated under negative pressure) 
to "funicular" (liquid is unsaturated and contiguous) , 
and finally to "pendular" (liquid is isolated) 35]. This 
drying process is mainly driven by the transportation of 
liquid water and vapor, and the phase change between 
them. In this paper, we adopt the water potential [36j as 
one of the key field variables. 

For the simplicity, we assume that the interactions be- 
tween the crack network and the water content field, the 
influences of the water pressure on the stress field, and 
the influences of the stress field on the water potential, 
indicated by dashed lines in Fig. [TJ are negligible. Our 
numerical simulations, however, exhibit good realizations 
of the experimental results as shown in the next sec- 
tion. The validity of these assumptions are discussed 
in Sec. IV. 

The stress field and its boundary are represented by 
the network of springs. And the water content field and 
its evolution are represented by the equations for the wa- 
ter potential. The details of the model are stated in the 
following subsections, i.e., the stress field, the water con- 
tent field using the water potential, and the interaction 
between them and their evolution. 



A. Stress field 

We construct a cubic lattice with nearest-neighbor 
(n.n.) and next-nearest-neighbor (n.n.n.) interactions in 
order to calculate the stress field [22| ■ Each interaction is 
represented by a Hookean spring. The lattice constant is 
clq and the size of the system is L x x L y x L z . The lattice 
point is expressed as a = (X, Y, Z), where < X < N x , 
< Y < N y , < Z < N z , L x = N x a , L y = N y a , 
L z = N z aQ. The z axis is chosen as the depth, the top 
(bottom) is set to z — (z — L z ). 

Let the spring constants for the nearest-neighbor and 
next-nearest-neighbor interactions be k\ and fc2, respec- 
tively. Then we obtain the spring network of the elastic 
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tensor Ciju with the components [2^ 



fci + 2fc 2 



a 



C, 



c • — I ■ 



«0 



(1) 

(2) 



We assume that the water content field changes very 
slowly and the stress field relaxes much faster than the 
water content field. So the stress field is assumed to be 
balanced at each moment under conditions induced by 
the water content field and the boundary condition at 
that time. We assume that each spring is broken at a 
certain critical force, i.e., if the force given on a spring 
exceeds a constant F c , the spring constant is set to be 
zero irreversibly. Using a position vector r a of the lat- 
tice point a, the force acting on the lattice point a is 
represented by 



(a a p - \r a ~ r/3 1) 



E * 



- _ r I ( a "/3~ l r a- r /3|): (3) 

r a ip 1 



A' 



/3£n.n.n. 

a/3 _ J fci (connected, i — 1,2) 
(broken) 



(4) 



where a a /3 is the natural length of the spring which con- 
nects lattice points a and (3. The equilibrium condition 
of the spring network, F a = 0, is equivalent to the con- 
dition that the elastic energy 



(5) 



where 



E — E a , 



E <* = E iKfia^-K-rpl) 2 

/SGn.n. 

+ E T^K"* 3 ( a <*P ~ \ r « ~ r p\f > (6) 

/3£n.n.n. 

is minimal under the given boundary condition and the 
given natural length of each spring. In our simulations, 
the free boundary condition is imposed. Because no fric- 
tion acts on the boundary of the system, the primary 
cracks do not appear. 



B. Water content field 

The evolution of the water content in unsaturated 
porous media consists of various processes. Here, we use 
the Campbell's desiccation model [37j that includes the 
transportation of liquid water and vapor and the evap- 
oration. By assuming that the temperature is constant 
and the speed of evaporation is slow, the heat of evap- 
oration is ignored. The evolution equation is described 



by one variable because both the liquid water content 
and the density of vapor are determined from the water 
potential as described below. 

The effect of the deformation of the spring network 
on the water content field is ignored and the evolution 
equation is solved on the undeformed lattice in the Carte- 
sian coordinates system because the change of the total 
volume is smaller than that of water volume (l3l | . The 
change of the boundaries due to the formation of cracks 
is also ignored. Note that there are many field variables 
and constants in this model. In order to distinguish them 
all the field variables are expressed by the symbols with 
an argument of position r (or z) and time t. 

The water potential is defined as the chemical poten- 
tial per unit mass of water in the mixtures, compared 
to that of pure, free water [H, [13, HI]. The potential is 
generally negative for unsaturated mixtures. Although 
this potential is the sum of several components, some 
components are negligible. The effect of the gravity are 
negligible because the size of a starch-water specimen in 
the experiments is on the order of 1 cm. We carried out 
the experiment [3^] to investigate the effect of gravity and 
confirmed that the columns normal to the free surface are 
formed even when the desiccated starch- water mixture in 
a container is laid on its side. We assume that there is no 
external pressure and the water in mixtures is pure, so 
we can ignore the overburden potential and the osmotic 
potential. Hence, the water potential is well represented 
by the matric potential ?/> m (r, t). The matric potential is 
defined as the amount of work per unit mass of water, 
required to transport an infinitesimal quantity of water 
from the mixture to a reference pool of water at the same 
elevation, pressure and temperature. The capillary water 
between particles and the water on the particle surface 
contribute to the matric potential. Below we use the 
matric potential as the water potential, i.e., the chemical 
potential per unit mass of water in mixtures. The contri- 
bution of the capillary water to the potential is expressed 



L cap 
PL 



(7) 



where P cap is the capillary pressure, expressed as 
—2o~/r cap , a is the surface tension of water and r cap is 
the inverse of the mean curvature of the capillary water 
in the mixture. 

From the coexistence condition of the liquid phase and 
the gas phase, the relative humidity h(r, t) is related to 
the matric potential ip m (r,t) by the Kelvin eq., 



h(r, t) = exp 



M w 
RT 



■0m(r,i) 



(8) 



where M w is the molecular weight of water, R is the gas 
constant and T is the absolute temperature. 

Let 0t(r, t) denote the volumetric water content. As- 
suming that 9l(y, t) = 9s in the saturated state {9s cor- 
responds to the porosity of the mixture), the volumetric 



4 



gas content 9c{r,t) is given as 

9 G {r,t) =9 s -e L (r,t). 



(9) 



During the drying process, 0j,(r,i) decreases as 9 G (r,t) 
increases. 9s is assumed to be constant because we ignore 
the feedback from the stress field (fracture and deforma- 
tion) to the water content field and the change of the 
pore volume due to the change of the water content. 

Expressing the fluxes of liquid water and vapor as 
Jt(r, t) and Jy(r, t) and the evaporation rate per unit 
volume as W eva (T,t), the transportation equation of liq- 
uid water and vapor is given by 



d_ 

Of 



(Mi(r,t)) = -V- J L {r,t)-W eva (r,t), (10) 



dt 



(pv(T,t)0 G (r,t)) = -V- J v (r,t) + W eva (r,t), (11) 



where pi, is the density of liquid water, pv( r , t) is the 
density of vapor. The incompressibility of water is as- 
sumed, i.e., pl is constant. 

In the unsaturated porous media, the flux of water is 
assumed to be proportional to the gradient of the water 
potential (Darcy's law) [11,13,11], 



J L {r,t) = -fci(r,t)Wm(r,i), 



(12) 



where fc£,(r,t) is the hydraulic conductivity. The prop- 
erties of starch-water mixtures in drying process, such 
as the dependences of the hydraulic conductivity &£,(r, t) 
and the water potential ?/> m (r, t) on the volumetric water 
content 9i{r,t) have not investigated in detail (40l. [4l|. 
Here, we suppose that the properties of starch- water mix- 
tures are similar to those of soil-water systems and em- 
ploy the relations confirmed by experiments in the soil 
mechanics 1351. 13711 . 



k L (r,t) = k L0 



9s 
h(r,t) 
9s 



2b+3 



where 



■0mO 

b 



-0.5d g 2 , 

1.0 x 10- 3 V>- 2 , 



0.2cr o 



(13) 
(14) 



(15) 
(16) 
(17) 



ipmo and fci,o are material constants. The units of ip m o 
and k^Q are J/kg and kg sec/m 3 , respectively. d g and a g 
are the geometric mean diameter and standard deviation 
of the particles, respectively (the unit is mm). The index 
b is a phenomenological parameter for soil-water systems. 
The relations (TIB")) and (TT4"|) are explained by some naive 
theories assuming the power-law distribution of pore size 
[37l ]. The dependences of fci(r,t) and |^ m (r,t)| on the 
diameter d g are indicated in Fig. [3J Further experimen- 
tal examinations may be required to confirm the validity 
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FIG. 2: Dependence of (a) kh and (b) \ip m \ on d g . A: d g = 



0.100 mm, B: d a 



0.050 mm, C: d a 



0.030 mm and D: 



for starch-water mixtures. We, however, expect that the 
important feature of this model is the functional shape 
of the two relations fci(t9i) and ip m (&L) and that the fol- 
lowing results are robust against tiny modifications. 
The flux of vapor is described by Fick's law, 



J v (r,t) = -D v e v e G (r,t)V Pv (r,t), 



(18) 



where Dy is the diffusion constant of vapor and ey is the 
tortuosity factor. The density of vapor is given by 



py(r,i) = p vo h(r,t), 



(19) 



where py is the saturation vapor density and the rela- 
tive humidity /i(r, t) is related to the matric potential by 
Eq. ©. 
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The evolution equation of tp m (r,t) is obtained from 
Eqs. 



= V • (fc ro (r, W m (r,*)), 



where 



at 



(20) 



(pz, - /9v(r,t))6>s / ip- 



fern ( r i t) 

fcy(r,t) 



-bip m (r,t) \ip m (r,t) 
M wPv (r,t)6 G (r,t) 
RT 

k L (r,t) + k v (r,t), 
D v e v G {r, t) Pv (r,t)M w 
RT 



(21) 
(22) 

(23) 



The boundary condition is that the evaporation occurs 
only at the top surface, 



kr 



dz 



-E, 



h(z = 0,t) - /i, 



-, (24) 



( Op 



aVf, 



a z 



= o, 

bottom 



(25) 
(26) 



where E p is the evaporation rate for the state that vapor 
is saturated on the top surface and h a is the atmospheric 
humidity and n is the normal vector at the side boundary. 
Note that the relative humidity at the top surface h(z = 
0, t) changes, so the evaporation rate is not constant [42j |. 



C. Combination of two fields 

The volume shrinkage occurs due to drying, which is 
the effect from the water content field to the stress field. 
This effect is incorporated into the model such that the 
natural length of each spring is an increasing function of 
the water content 0£,(r, t), i.e., the natural length a a p of 
the spring between the lattice points a and (3 is deter- 
mined as 



a a /3 



a.0a(3 <l + K 



6lo + &L/3 



\N X ' N y N z ' 



f a (n.n.) 

n.n.n.) 



(27) 
(28) 
(29) 



The parameters of the stress field, fci, &2, the expansion 
coefficient k, and the critical force for the breaking F c , are 
supposed to be constant for the simplicity. The changes 
in the stress field may affect the water content field (or 
the water potential) , however this effect is assumed to be 
ignored. 

We also assume that evaporation does not occur from 
the crack surfaces as well as the side boundaries. There- 
fore, in our model, all the field variables, such as the 



water content field 0£,(r, t) and ^> m (r, t) are the function 
of the depth z and time t. The argument of the field 
variables are changed from (r,t) to (z,t) hereafter. The 
drainage effect, i.e., the role of the evaporation through 
the crack will be discussed in Sec. IV. 



D. Time evolution 

The initial states in our simulations are prepared so 
that the pore is saturated, 



■Az,t = Q) = 6 s , 



(30) 



and randomly chosen 1% of the springs at the surface(z = 
0) are broken. 

We repeated the following procedures at each time 
step: 

1. Time steps by dt and the water potential ip m (z,t) 
is calculated from Eq. (f2"0"|) and the water content 
is obtained from Eq. 11140. 



2. The natural length of each spring is determined 
from the water content field 0L(z,t) by Eqs. (j2"T|) - 



3. The equilibrium configuration of the stress field is 
calculated so that the elastic energy ^ has the 
minimal value by the conjugate gradient method 
[43| | with a tolerance 10~ 4 . 

4. If the force given on a spring exceeds the critical 
value F c , it breaks, i.e., the corresponding spring 
constant changes to zero. If more than one spring 
have loads larger than F c , break one which has the 
maximum force. 

5. The procedures 3 and 4 are repeated until no spring 
is broken, and then return to the first procedure. 



III. NUMERICAL RESULTS 
A. Typical patterns 

The parameters used in the simulations are displayed 
in Table[l] The elastic and fracture parameters ki, ki and 
F c are chosen as it works. Below we show the results for 
the particle diameter d g = 0.050 mm except for Figs. [H 
and [T31 The typical results at t = 100 h which corre- 
sponds to the middle of the stage 3 are shown in Fig. [3] 
The water content decreases due to drying and the region 
where the water content field changes sharply exists, i.e., 
the drying front exists near z ~ 0.4 cm. The non-uniform 
contraction near the drying front consequently increases 
the stress in this region and fractures occur. In the front 
of cracks, the stress intensifies and the tensile springs (de- 
noted by red or dark gray in Fig. [3}}) accumulate. Com- 
pressive springs (blue or light gray) accumulate between 




FIG. 3: (Color online). Typical results of the numerical sim- 
ulations. The unit of axis is cm. (a) The water content field 
9l at t = 100 h. (b)The stress field at t = 100 h. The springs 
of which force is more than 0.4F C are plotted in red (tensile) 
and blue (compressive), or for grayscale image, in dark gray 
(tensile) and light gray (compressive). In the front of cracks, 
the stress intensifies and the tensile springs accumulate. Com- 
pressive springs accumulate between the crack tips. (c)The 
crack pattern. The springs broken in 90 < t < 100 h are 
shown. 
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system size L x x L y x L z 


2.5 x 2.5 x 2.5 cm 3 


N x x N y x N z 


50 x 50 x 50 


density of liquid water pl 


1.0 x 10 3 kg m" 3 


density of saturation vapor pvo 


1.7 x 10~ 2 kg m~ 3 


saturation water content 8s 


0.5 


mass of a mol of water M w 


18 g mol" 1 


gas constant R 


8.3 J mol" 1 K" 1 


temperature T 


293 K 


RT/M W 


1.35 x 10 5 J kg" 1 


diffusion constant of vapor Dv 


2.5 x 10" 5 m 2 s _1 


tortuosity factor ev 


0.66 


evaporation rate E v 


0.01 g cm -2 h _1 


atmospheric humidity h a 


0.5 


geometric mean particle diameter d g 


0.01 ~ 0.10 mm 


geometric standard deviation 


1.0 


of a particle diameter, a g 




spring constants k\ao, feao 


1.0, 0.5 


critical force for the breaking F c 


0.003 


expansion coefficient k 


0.2 


time step dt 


0.1 h 



TABLE I: Parameters of the water content field [§3] and the 
stress field. 



the crack tips. Figure [3{c) represents the springs broken 
between t = 90 h and t = 100 h. The region where frac- 
tures occur actively forms a polygonal pattern and the 
trace of such region develops into the prismatic structure 
of the cracks. These results well capture the features ob- 
served in the real experiments [HI, HH, [H, [l?], EU ■ Next, 
we focus on the spatio-temporal behavior of the key el- 
ements, i.e., the water content field, the stress field and 
the crack pattern. 



B. Water content field 

The time evolution of the normalized total mass of liq- 
uid water J 8l{z, t)/0sdz is displayed in Fig. [4] The two 
different stages are observed as in the experiments [l3| . 
i.e., constant slope stage and slower drying stage. The 
former and the latter correspond to the stage 1 and 3, 
respectively, and the stage 2 connects them. The time 
when the stage 1 switches to the stage 3 depends on the 
particle diameter d g . From the Kelvin eq. ([5]) and the 
relation between the water content and the water poten- 
tial, Eq. (jT5J) (shown in Fig. as the particle diameter 
is large, the large amount of water evaporates. 

The depth dependences of the several functions related 
to the water content field are depicted in Figs. l5l6t71 and[51 
In the stage 1 (0 < t < 50 h), the volumetric water con- 
tent 6l decreases uniformly and the flux of liquid water 
is dominant to the transportation of water. In the stage 
3 (t > 80 h), the drying front appears and propagates 
inward. The dominant process of the transportation of 
water differs between above and below the front, i.e., the 
flux of liquid water Jl dominates below the drying front 
while the flux of vapor Jy replaces it (Figs. [S]and[n]). 
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50 100 150 200 250 

Time (h) 



FIG. 4: Evolution of normalized total mass of liquid 
water/ 6 L (z,t)/6 s dz. A: d g = 0.100 mm, B: d g = 0.050 mm, 
C: d g = 0.030 mm and D: d g = 0.010 mm. 
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FIG. 6: Flux of liquid water Jl and flux of vapor Jv vs depth 
z. L : J L (t = 50h), Li: J L {t = lOOh), Vi: J v (t = lOOh), L 2 : 
J L (t = 150h) and V 2 : Jv(t = 150h). J v (t = 50h) is on the 
order of lO" 11 Kg m~ 2 s -1 . 




Above the drying front, the absolute value of the wa- 
ter potential \ip m (z,t)\ increases sharply to the order of 
RT/M W (~ 10 5 J/kg) and the relative humidity decreases 
below significantly 1 (Figs. [7] and [5]). The contribution 
of the capillary water between micrometer-size particles 
to the water potential is on the order of 10 2 J/kg. So, 
the effects of the electrical and van der Waals force in 
the vicinity of the solid-liquid interface are supposed to 
be important (38j . 



C. Stress field 

Figure [9] displays the water content field, the normal- 
ized elastic energy density and the stress field on the 
vertical section of the system at t = 120 h. The en- 
ergy density is normalized by the maximum value in the 
vertical section. There are the drying front where the 
water content sharply changes and the region where the 
tip of cracks exists and the region where tensile springs 
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(denoted by red or dark gray in Fig. Etc)) accumulate. 
Compressive springs (blue or light gray) accumulate be- 
tween the crack tips and stress does not concentrate in 
the region far from the drying front. The trace of the frac- 
tured region forms the columnar crack pattern observed 
finally. There is a crack of which tip does not follow to 
the drying front [x ~ 2.0 cm). It is an elementary process 
of columnar merger. 



D. Cracks 

Cracks are formed when the uniformity of the water 
content distribution is lost and the tips of cracks prop- 
agate inward. The time evolution of X)r<t^( 2 ' r ) are 
shown in Fig. IIOI where l(z, r) represents the number of 
springs broken at (z, r). Y] T<t l(z, r) corresponds to the 
cumulative crack length at (z,t). The area surrounded 
by the two curves at different times corresponds to the 
number of springs broken between those times. From the 
fact that the difference between the distribution of broken 
springs for two adjacent time has finite support, the ac- 
tive cracking area is localized in z-direction. The relation 
between the position of active cracking area and the dry- 
ing front can be checked as follows. The space-time plots 
of d0i,(z,t)/dz, ^y£ a and l(z,t), are shown as gray 
scale images in Fig. [n] The energy density is normalized 
by the maximum value in all region. It is confirmed that 
the active cracking area moves inward following the dry- 
ing front. However, the number of the broken springs is 
not the simple function of the water content 6l{z,{) or 
its derivative. 

The typical patterns of cracks formed up to the final 
state it — 200 h) are shown in Fig. [T2] The polygonal 
columnar structure can be seen despite of the anisotropy 
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FIG. 9: (Color online). The unit of axis is cm. (a)The 
water content field, (b)the normalized elastic energy density, 
and (c)the stress field are displayed in 1.0 < y < 1.25 cm at 
t — 120 h. The springs of which force is more than 0.3F C 
are plotted in red (tensile) and blue (compressive), or for 
grayscale image, in dark gray (tensile) and light gray (com- 
pressive). The broken springs are plotted in green, or for 
grayscale image, in light gray. In the drying front where the 
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FIG. 10: Time evolution of the cumulative number of the 
broken springs at each horizontal cross sections. The results 
obtained from seven different initial conditions are averaged 
at a given depth z. 



of the cubic lattice. In ref. [la ], two mechanism for chang- 
ing scale are suggested, i.e., merging and column initia- 
tion from a vertex. In Fig. [T2l both process are observed. 
These features well capture the exp eriments of starch- 
water mixtures [H, [lj| [H, [13, O- Now, we analyze 
the depth dependence of the final crack pattern. Fig- 
ure [TS] exhibits the depth dependence of the cumulative 
crack length up to the final state with double logarithmic 
plot to compare with the experimental data. In contrast 
to the real experiments, we can easily change the phys- 
ical parameters. If we change the particle diameter d gi 
it affects the VmO through phenomenological Eq. ([T5J) . 
The size of columns decreases as the particle diameter 
increases. Though it is difficult to state definitely, there 
is a regime in which the number of broken springs de- 
creases with a power law with exponent between —0.5 
and —2. 



IV. DISCUSSIONS 

In this section, we note several points such as the trans- 
lation of the evolution equation from the matric potential 
-0 m (z, t) to the water content 9l(z, t), the drainage effect 
of cracks, the geological columnar joints, the feedback 
from the stress field to the water content field, the char- 
acterization of the crack patterns and the properties of 
the granules water mixtures. 

First, we discuss the drying front which is obtained 
from the evolution equation of the water potential 
(Eq. ([2l7)) V Changing the variable from -0 m (z, i) to 
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FIG. 11: Space-time plot of (a) d0L/dz, (b) the normal- 
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springs. These fields are averaged over the cross section at a 
given z. The line indicates the depth where the flux of vapor 
is equal to that of liquid water (see Fig. [6} , which corresponds 
to the position of the drying front. 
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6> L (z,i) by Eq. (H}, we get 
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FIG. 13: The cumulative number of the broken springs at 
each cross section at the final state (t = 200 h) for various 
particle diameters, A: d g = 0.100 mm, B: d g = 0.050 mm, C: 
d g — 0.030 mm and D: d g — 0.010 mm. The results obtained 
from seven different initial conditions are averaged for each 
particle diameter. 



where 



D L (z,t) 



Pv(z,t) / b^ m (z,t)M w 9 G (z,t) 
PL \ RT9 L (z,t) 

bip m {z, t)(k L (z, t) + k v (z, t)) 



e L {z,t) 



(32) 
(33) 



(31) 



g{9L) is of order 10~ 4 due to pv <C Pl\ therefore, if we ig- 
nore g{0i)t the time evolution of 0l can be considered to 
obey the nonlinear diffusion equation with variable diffu- 
sion coefficient Dl(z, t). In our simulations the diffusion 
coefficient Di(z,t) has a minimum near the drying front 
(Fig. [14]) as Goehring et al. [18| pointed in their ex- 
periments. The nonlincarity of the diffusion coefficient 
is crucial for the water content distribution which has a 
concave regime. 

Second, we discuss the assumption that cracks do not 
affect the water content field. Assuming the crack aper- 
ture width l c r~j 10 _1 mm and the propagation speed of 
the drying front v w ~ 10~ 5 mm/sec and the volumetric 
gas content near the drying front Og ~ 0.4, the length 
of the vapor diffusion in t w — l c /v w is on the order of 
y/DyevOcTw ~ 10 2 mm. Because this is much larger 
than the crack aperture width, the density of vapor in a 
crack aperture is determined quasi-statically by the wa- 
ter potential on the crack surface. Therefore, the effect 
of cracks on the dynamics of the water content field as 
the boundary condition can be negligible above the dry- 



11 




Depth (cm) 

FIG. 14: Diffusion coefficient D L defined by Eq. (J33J. 



ing front because the flux of vapor is dominant in this 
region. The flux of liquid water may be affected by the 
cracks in the region where the two fluxes are comparable. 
However we infer that the flux of water does not change 
significantly because the region is limited in the vicinities 
of the crack tips. Hence the assumption is valid except 
for the small effect near the drying front. 

For studying the effect of cracks on the water content 
field, we also perform numerical simulations based on the 
different assumptions that the transportation of liquid 
water obeys the linear diffusion equation and the diffu- 
sion constant in cracks is larger than in medium. The 
result, however, is that the broken springs do not form 
planes to become crack surfaces but form some clusters 
and no stable prismatic pattern is observed. Recently, 
S0renssen et al. [25| report a similar fracture simulation 
under the condition that the quantity which determines 
the local volume reduction diffuses to the boundary of 
the solid and find that the zone where the fracture occurs 
actively propagates and the fracture patterns resemble a 
forest of trees, i.e., there is many branches of cracks. In 
order to obtain the prismatic fracture pattern under the 
condition that the water content field does not affect the 
stress field except for determining the volume shrinkage 
rate, the mechanism that the external field forms the 
front by itself without cracks would be necessary. It is 
the nonlinear diffusion of water in our simulations. 

Third, let us consider the another columnar joint for- 
mation phenomenon, i.e., the geological columnar joint 
formation driven by lava cooling process. It is accepted 
that the cooling process of lava is the combination of 
the linear diffusion of heat and the convective heat loss 
through cracks @, 0, S HE H3. On the other hand, 
in our spring network model we obtained the prismatic 
structure of cracks in the condition that cracks do not 



affect the water content field and the water pressure do 
not affect the stress field. There is the difference in the 
external field, i.e^ the thermal field above the front is 
nearly constant [26[ and the water content field below the 
front is nearly constant [l?], EH ■ So it is a future work 
to study numerically the cooling joints with taking ac- 
count of the convective heat loss through cracks and the 
external stresses i.e., the horizontal tectonic stress, the 
vertical overburden stress and the local pore pressure. 

Fourth, we discuss the feedback from the stress field 
to the water content field. We ignore this feedback, i.e., 
the water potential is assumed not to be affected by the 
change of the stress field. In our simulations, as described 
in Sec. Ill B, the contribution of the capillary water pres- 
sure to the water potential becomes much smaller than 
the absolute value when the water content decreases. 
Therefore, we infer that, when the drying front appears, 
the deformation does not change significantly the water 
potential through the water pressure. This inference is 
consistent with our assumption. However it is possible 
that the deformation may affect the properties of water 
on the particle surfaces and the shapes of the water- vapor 
interfaces. It is a future subject to investigate the con- 
tribution of these effects to the water potential. 

Fifth, we discuss the characterization of the crack pat- 
terns. The fracture occurs near the propagating front 
of the water content field. The trace of the fractured 
region forms the columnar pattern with polygonal cross 
sections. The size of columns increases with depth. In 
this paper, as a key quantity, we focus on the crack length 
(i.e., the number of broken springs) which is supposed to 
be proportional to the surface energy in connection with 
the Griffith theory [44[. Other characterizations, such 
as the mean polygon area and number of polygon sides, 
seems to be worth analyzing as a future work. However, 
they should require larger system size and careful defini- 
tion of the quantities. If the pattern at the cross section 
is assumed to be isotropic, the characteristic length of the 
pattern is inversely proportional to the crack length, so 
it is possible to compare the depth dependence of the cell 
area at the cross section with the experimental results. 
We suggest that the drying granules-water mixtures of 
larger particles forms the smaller size of columns, the 
dependence of the size of columns on the particle diam- 
eter is a future subject. 

Finally, we discuss the properties of the granules wa- 
ter mixtures. We obtain the prismatic pattern of cracks 
by employing the relations for general soils in our simula- 
tions. For the present, it is known that only the combina- 
tion of starch and water forms the prismatic fracture pat- 
tern among drying mixtures of various combinations of 
granules and liquids. For understanding why dried mix- 
tures of other granules and water do not form the colum- 
nar structure, it is necessary to measure experimentally 
various material properties such as the expansion coeffi- 
cient and the dependences of the water potential and the 
hydraulic conductivity on the water concentration. 
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V. CONCLUSIONS 

We study the drying process and the columnar frac- 
ture process in granules-water mixtures. The mixture is 
assumed to be the elastic porous medium and each pro- 
cess is represented by a phenomenological model with the 
water potential and a spring network, respectively. We 
assume that the cracks do not affect the water content 
field as a boundary condition and ignore the feedback 
from the stress field to the water content field. In the 
numerical simulations, we find that the water content 
distribution has a propagating front and the pattern of 



cracks exhibits a prismatic structure. 
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